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EFFECTS  OF  CAVITATION  ON  UNDERWATER 
SHOCK  LOADINGS  -  Part  I 

Abstract 

Reported  here  are  analytic  formulations,  together  with 
one-dimensional  results,  in  an  investigation  of  the  title 
subject.   It  is  shown  that  either  displacement  or  a  displacement 
potential  may  be  used  as  the  basic  dependent  variable  for  a 
finite  element  analysis.   Artificial  damping  is  found  to  be 
needed  to  suppress  spurious  oscillations  (a  numerical  phenom- 
enon) near  cavity  boundaries.   Adequacy  of  the  method  is  demon- 
strated by  comparison  with  published  results  of  Bleich  and 
Sandler.   Some  results  are  given  for  effects  of  cavitation  on 
the  performance  of  resilient  attenuators. 


EFFECTS  OF  CAVITATION  ON  UNDERWATER 
SHOCK  LOADINGS  -  Part  I 

1.   Introduction 

1.1  Earlier  Work.  Motivation  for  the  present  investigation  origin 
ated  with  the  "shock  shield"  proposed  by  Geers  in  Ref.  1.   The 
shield  is  a  gas-filled  cushion  (GFC)  to  be  fitted  to  the  exterior 
of  a  submarine  hull.   If  the  cushion  has  sufficient  thickness 

it  can  greatly  reduce  the  magnitude  of  underwater  shock  loads 
transmitted  to  the  hull.   Related  concepts  involving  the  applica- 
tion of  resilient  elastic  layers  (REL)  are  treated  by  Geers  in 
Ref.  2. 

1.2  Cavitation  Effects.   The  two-dimensional  analyses  of  Refs.  1 
and  2  neglect  possible  effects  of  cavitation.   It  is  well-known 
(see  Refs.  3  and  4),  however,  that  a  highly  compliant  submerged 
object  will  produce  a  negative  pressure  scattered  wave  in  response 
to  an  incident  shock  wave.   If  the  shock  pressure  is  much  greater 
than  the  hydrostatic  pressure,  cavitation  will  be  induced  in  the 
fluid.   Such  cavitation  may  significantly  increase  the  shock  load- 
ing on  the  submerged  body.   It  is  thus  evident  that  an  adequate 
investigation  of  the  effectiveness  of  resilient  attenuators 
requires  evaluation  of  effects  of  cavitation  on  performance. 

1.3  Plan  for  Investigation.   The  present  investigation  is 
divided  into  two  phases.   The  first  phase  is  the  subject  of  this 
report.   It  involves  consideration  of  representative  one- 
dimensional  problems  for  the  purpose  of  determining  the  relative 
merits  of  alternate  choices  for:  dependent  variables,  time 


integration  algorithms,  and  spatial  and  temporal  discretization. 
The  one-dimensional  context  allows  rapid  and  inexpensive  compu- 
tations and  allows  comparison  of  results  with  those  reported  by 
others.   An  account  of  this  phase  is  given  in  the  remaining 
sections  of  this  report. 

The  second  phase  of  the  investigation  consists  of  extensions 
to  problems  in  three  dimensions.   Reasonable  limitations  on  com- 
puter core  capacity  and  processing  time  require  that  the  problems 
be  axisymmetric  and,  thus,  mathematically  two-dimensional. 

1.4   Fluid  Model .   It  is  known  that  fluids  do  have  some  capacity 
for  sustaining  negative  pressure  (tension) .   Some  data  are  given 
in  Ref.  5.   The  influence  of  dissolved  gas  on  the  development 
of  cavitation  is  considered  in  Ref.  6.   For  the  purpose  of  the 
present  investigation  it  is  advantageous,  and  presumably  con- 
servative, to  assume  that  the  transition  from  the  normal  to  the 
cavitated  state  takes  place  without  delay  when  the  absolute 
pressure  reaches  zero. 

In  the  initial  stages  of  this  investigation  the  fluid  was 
treated  as  bilinear  with  a  greatly  reduced  bulk  modulus  in  the 
negative  pressure  region.   Subsequent  developments  disclosed 
that  the  expected  advantages  of  the  bilinear  model  were  not 
achieved  and  the  bulk  modulus  was  henceforth  assumed  to  be  zero 
in  the  cavitated  region. 


2.   Choice  of  Dependent  Variable 

2.1  Failure  of  the  Pressure  Formulation.   At  the  outset,  this 
investigator  expected  that  a  formulation  of  the  governing 
equations  using  fluid  pressure  p  as  the  basic  dependent  variable 
would  be  advantageous.   This  expectation  was  based  on  previous 
successful  finite  element  applications  to  propagation  problems 
(e.g.,  see  Ref s .  7-9).   Prior  applications  did  not  involve  cavi- 
tation, but  the  bilinear  fluid  model  was  expected  to  handle 
successfully  cavitation  effects. 

At  an  early  stage  of  the  investigation,  duplication  of  the 
results  of  the  example  problem  of  Ref.  4  was  attempted.   These 
trials  gave  solutions  which  correctly  tracked  the  growth  of  the 
cavitated  region,  but  failed  to  show  its  subsequent  contraction 
and  collapse.   Efforts  to  discover  the  reason  for  the  failure  of 
the  pressure  formulation  led  to  a  simple  test  problem  which 
determines  whether  a  proposed  formulation  can  correctly  track 
the  contraction  of  a  cavitated  region.   The  problem  is  defined 
in  the  following  section. 

2.2  "Water-Hammer"  Problem.   The  rapid  pressure  rise  which 
accompanies  the  sudden  interruption  of  water  flow  in  a  closed 
conduit  is  known  as  water_hammer .   We  here  consider  a  flow  in  a 
zero  pressure  (or  a  small  negative  pressure,  if  the  bilinear 
fluid  model  is  used)  cavitated  region  with  positive  dilatation 

e  .   Thus  we  have,  for  a  semi- infinite  region  x  >  0,  the  initial 

values : 

p(x,0)  =  0,   (pressure) 

e(x,0)  =  e0  >  0,   (dilatation) 
u(x,0)  =  -v  <  0.   (velocity) 


The  boundary  condition  at  x  =  0  is  u(0,t)  =  0.   The  exact 
solution  to  this  problem  is  especially  simple.   A  shock  front 
propagates  with  constant  speed  etc,  beginning  at  the  closed  end, 
and  the  fluid  behind  the  front  is  at  rest  with  uniform  pressure 
p,  =  apev  .   Meanings  of  the  symbols  introduced  are: 

p  =  fluid  density, 
c  =  acoustic  velocity. 
Factor  a  is  given  by 

a  =  [1  +  (ce  /2v  )2]h    -  (ce  /2v  ) . 
u     *•   o    O  0    0 

In  the  region  ahead  of  the  shock  front  (x  >  act)  the  variables 
p,  u,  and  e  maintain  their  initial  values. 

2.5   Governing  Equations.   Details  concerning  the  governing 
equations  are  given  in  Appendix  A.   Equations  are  stated  in 
forms  suitable  for  any  number  (i.e.,  one,  two,  or  three)  of 
spatial  dimensions.   Four  separate  formulations  are  developed 
with  the  principal  dependent  variable  being  particle  displacement, 
fluid  pressure,  velocity  potential,  and  displacement  potential 
for  the  respective  cases.   The  capability  of  each  to  deal  with 
the  one-dimensional  water-hammer  problem  is  discussed  separately 
below. 

2.4   Displacement  Formulation.   In  this  case  it  is  convenient  to 
replace  the  displacement  vector  <5_  by  r  =  p5_  and,  instead  of  the 
dilatation  e,  a  density  -  weighted  condensation  s  =  -pe  is  intro- 
duced.  The  working  equations  (Eqs.  All,  A12,  and  A10)  consist 
of  an  equation  expressing  r  in  terms  of  the  gradient  of  p,  one 


giving  s  as  the  negative  of  the  divergence  of  r,  and  the  pair  of 
algebraic  relations  for  calculating  p  from  s  (the  bilinear  con- 
stitutive law).   It  is  readily  seen  that  the  initial  conditions 
of  the  water-hammer  problem  allow  determination  of  the  initial 
values  of  r  and  r.   The  first  may  be  deduced  through  spatial 
integration  of  the  constant  initial  dilatation  e  and  the  second 
is  known  from  the  given  initial  velocity  -v  .   Given  these 

a  j  q 

required  initial  values  and  the  boundary  condition  at  x  =  0 , 
the  working  equations  suffice  to  solve  the  problem. 

The  displacement  vector  is  continuous,  but  both  velocity 
and  pressure  are  discontinuous  at  the  shock  front.   The  pressure 
discontinuity  is  not  representable  by  the  shape  functions  of  the 
finite  element  method.   Accordingly  there  must  be  a  finite 
length  interval  over  which  both  the  pressure  rise  and  particle 
deceleration  take  place.   The  necessity  for  this  compromise  must 
be  considered  a  defect  (but  not  a  disqualifying  one)  of  the  dis- 
placement formulation.   A  further  disadvantage,  not  shared  by 
any  of  the  other  formulations  is  that  the  principal  dependent 
variable  is  a  vector,  not  a  scalar.   In  the  axisymmetric  appli- 
cations planned  this  will  double  the  order  and  bandwidth  of  the 
stiffness  matrix:,  resulting  in  a  great  increase  in  requirements 
for  computer  storage  and  processing  time. 

2,5   Pressure  Formulation.   In  Appendix  A  the  pressure  formulation 
is  stated  in  terms  of  a  second  order  equation  (A  13)  expressing 
s  as  the  Laplacian  of  the  dynamic  component  of  pressure,  together 
with  the  bilinear  constitutive  relation  (A  10).   The  solution 
technique  involves  stepwise  time  integration  for  s,  alternating 


with  use  of  the  constitutive  relation  to  find  new  values  of  p. 
In  view  of  the  role  played  by  s,  it  is  slightly  misleading  to 
call  this  the  pressure  formulation.   Indeed,  in  the  absence  of 
cavitation  it  is  advantageous  to  use  the  constitutive  relation 
to  eliminate  s  in  favor  of  p  alone.   In  cavitated  regions,  how- 
ever, retaining  s  makes  the  strategem  of  a  nonzero  bulk  modulus 
unnecessary.   Despite  this  advantage,  both  versions  fail. 

The  reason  this  technique  fails  when  applied  to  the  water- 
hammer  problem  is  readily  apparent.   The  variable  s  (and  also  p) 
has  a  step  discontinuity  at  the  shock  front.   (Correspondingly, 
the  second  derivative  s  has  a  dipole  singularity.)  We  know 
that  the  height  of  the  step  depends  on  the  initial  (negative) 
value  of  s  and  on  the  fluid  velocity.   The  velocity,  however, 
is  neither  explicitly  nor  implicitly  represented  in  the  equations 
There  seems  to  be  no  further  reason  to  consider  the  pressure 
formulation  for  cavitation  problems. 

The  bilinear  fluid  model  with  a  nonvanishing  bulk  modulus 
in  the  tension  region  was  initially  introduced  with  the  expec- 
tation that  the  pressure  formulation  would  work.   Since  the 
sole  reason  for  including  this  complication  has  vanished,  we 
shall  henceforth  exclude  negative  pressures  and  correspondingly 
choose  3  =  0  in  (A2) ,  (A4) ,  and  (A10) . 

2.6   Velocity  Potential  Formulation.   Using  4>  to  represent  the 
velocity  potential,  the  governing  equations  express  <j>  as  the 
negative  of  the  dynamic  pressure  (A14) ,  s  as  the  negative  of  the 
Laplacian  of  $  (A15) ,  and  p  in  terms  of  s  (A10) .    In  this 


instance  the  velocity  potential  is  continuous,  but  its  first 
derivatives  are  not.   It  is  possible  to  obtain  a  moderately 
satisfactory  solution  from  the  discretized  equations. 

2.7  Displacement  Potential  Formulation.   Using  \p    to  represent 
the  displacement  potential,  the  governing  equations  express  ip 

as  the  negative  of  the  dynamic  pressure  ( A 1 6 ) ,  s  as  the  negative 
of  the  Laplacian  of  ip    (A17)  ,  and  p  in  terms  of  s  (A10)  .   All  of 
the  needed  variables  appear,  explicitly  or  implicitly,  in  this 
formulation.   Moreover,  ip  and  its  first  derivatives  are  contin- 
uous.  The  step  discontinuities  in  p,  s,  and  u  are  manifested 
as  discontinuities  in  the  second  derivatives  of  ty . 

2.8  Formulation  Selection.   Among  the  formulations  considered, 
the  pressure-based  one  is  discarded  as  unworkable  and  the  velo- 
city potential  is  rejected  as  inferior  to  the  displacement  poten- 
tial on  the  basis  of  continuity.   In  the  remainder  of  this  report, 
discussion  of  application  details  will  be  confined  to  the  \p 
formulation  because  it  is  novel.   The  displacement  formulation 
was  developed  in  parallel  and  also  tested  on  the  Bleich-Sandler 
example  (Ref .  4) . 

2.9  Discretized  Equations  and  Solution  Process.   The  process  of 
forming  discretized  finite  element  equations  from  the  correspond- 
ing partial  differential  equation  is  well  known  (e.g.,  see  Ref.  7) 
and  will  not  be  repeated  here.   We  note  that,  based  on  prior 
experience  with  wave  propagation  problems,  linear  shape  functions 
were  chosen.   It  was  found  that  a  lumped  "mass"  matrix  was  easier 
to  use  and  gave  better  performance  than  its  consistent  counterpart 


Details  concerning  the  formulation  of  initial  conditions 
and  the  radiation  boundary  condition  are  given  in  Appendix  B. 
The  considerations  used  for  selecting  a  time  integration  algor- 
ithm and  the  means  for  introducing  needed  artificial  damping 
are  discussed  in  Appendix  C. 

3.   The  Bleich-Sandler  Example 
3.1   Statement  of  the  Problem.   In  Ref.  4  Bleich  and  Sandler, 
using  a  bilinear  fluid  model,  study  cavitation  phenomena  during 
one-dimensional  wave  propagation.   They  use  the  method  of 
characteristics  and  introduce  additional  relations  to  connect 
state  variables  on  opposite  sides  of  a  shock  front. 

The  numerical  example  they  give  concerns  "the  response  of  a 
horizontal  layer  of  mass  on  the  surface  of  a  half-space  of  fluid, 
Fig.  1.   A  plane  pressure  wave  with  a  sudden  rise  and  an  expon- 
ential decay  moves  toward  the  surface,  reaching  the  mass  at  the 


Surface  mass    W 

r 


Atmospheric 
pressure 


■  ace  mass    w  •,  pressure 


Fluid 

half-  space 

x  >0 

Pressure  History 

Figure  1.   Particulars  of  Bleich-Sandler  example  (from  Ref.  4). 

time  t  =  0.  The  system  is  subject  to  gravity  and  atmospheric 
pressure,  all  particles  being  at  rest  prior  to  arrival  of  the 
shock.   The  analysis  is  based  on  the  degenerate  model  with  S  =  0." 


10 


For  calculation  and  presentation  of  results  Bleich  and  Sandler 
use  the  time  constant  of  pressure  wave  decay  (~  1  ms)  as  the 
time  unit.   The  acoustic  velocity  c  is  given  the  convenient 
value  unity  by  choosing  unit  length  to  be  the  distance  traveled 
by  the  pressure  wave  in  unit  time  (-  56  in.). 

3.2   Comparisons  with  Bleich-Sandler  Results.   Bleich  and 
Sandler  present  two  figures  summarizing  their  solution.   First 
of  these  traces  the  time  history  of  the  cavitated  region  in  the 
x  -  t  plane.   Their  figure  is  reproduced  in  Figs.  2  and  3  below 
with  superposed  points  obtained  by  present  analyses.   For  the 
finite  element  analyses  the  discretized  region  extends  from 
x  =  0  to  a  radiation  boundary  at  x  =  4.   Results  shown  in  Fig.  2 


t,  non-dimensional  time 


Case? 


Figure  2. 


Bleich-Sandler  example:  time-history  of  the 
cavitated  region.   Discrete  points  found  by 
finite  element  method  usirg  the  displacement 
formulation. 
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were  obtained  from  the  displacement  formulation  using  a  node 
spacing  (element  length)  Ax  =  .04  and  a  time  step  h  =  .01.   In 
the  absence  of  damping  it  was  found  that  moving  "islands"  of 
positive  pressure  appeared  within  the  cavitated  region.   This 
behavior  is  henceforth  called  "frothing."   It  was  found  that 
damping  with  n  =  .16  was  sufficient  to  supress  frothing  and 
produce  a  smooth  variation  of  the  condensation  s  within  the 
cavitated  region. 

In  Fig.  3  the  superposed  points  were  obtained  using  the 
displacement  potential  formulation.   For  these  results:  Ax  =  .01 


t,   non-dimensional  time 


Case  7 


Figure  3.   Bleich-Sandler  example:  time-history  of  the  cavitated 
region.    Discrete  points  found  by  finite  element 
method  using  the  displacement  potential  formulation. 
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and  h  =  .0025.   To  suppress  frothing  a  value  ri  =  .0025  was 
found  to  be  sufficient. 

A  second  Bleich-Sandler  figure  shows  the  time  history  of 
the  velocity  of  the  surface  mass  layer.   This  is  reproduced  as 
Figs.  '4  and  5.   Points  obtained  from  present  analysis  using  the 
displacement  formulation  are  superposed  on  Fig.  4  and  points 


) 

12  • 


I1-- 


--ts,  arrival  of  secondary  shock  al  the  surface 
--tw,  closure  of  cavited  region 


with  cavitation 


1000  u. 


01  0  2  0.3  0.4  0  5 


Figure  4.   Bleich-Sandler  example:  nondimensional  upward  velocity 
of  surface  mass.   Discrete  points  found  by  finite 
element  method  using  the  displacement  formulation. 


from  the  displacement  potential  formulation  on  Fig.  5.   The 
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t_,  arrival  of  secondary  sriock  at  the  surface 


--tw,  closure  of  cavited  region 


1000  u. 


0.1  0.2  0  3  0.4  0  5 


Figure  5.   Bleich-Sandler  example:  nondimensional  upward 

velocity  of  surface  mass.   Discrete  points  found 
by  finite  element  method  using  the  displacement 
potential  formulation. 


parameters  used  were  those  given  for  Figs.  2  and  3,  respectively 
It  is  believed  that  the  agreement  demonstrated  in  Figs.  2-5 
substantiates  the  adequacy  of  displacement  and  displacement 
potential  formulations  for  one -dimensional  analyses. 
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4.   Effect  of  Cavitation  on  Resilient  Attenuator  Performance 

We  consider  here  the  effect  of  cavitation  on  the  perform- 
ance of  two  selected  attenuators.   One  is  a  GFC  with  an  initial 
thickness  of  10  in.  under  hydrostatic  pressure.   The  second  is 
an  REL  with  modulus  C-.  =  1000  psi.  and  an  initial  thickness 
L   =10  in.  when  loaded  only  by  atmospheric  pressure  .   The 
shock  loading  consists  of  a  pressure  step  of  amount  p   followed 
by  exponential  decay  with  time  constant  25  ms .   The  hydrostatic 
pressure  is  p,  .   Results  are  given  as  the  quotient  of  the  maxi- 
mum dynamic  pressure  increment  Ap    by  the  peak  shock  pres- 
sure p  . 
rs 

If  cavitation  effects  are  ignored  the  response  may  be  found 
by  integrating  numerically  a  first  order  ordinary  differential 
equation.   Details  are  omitted  here. 

The  response  considering  cavitation  has  been  determined 
using  finite  element  modelling  based  on  the  displacement  poten- 
tial formulation.   Results  are  summarized  in  Table  1. 

Table  1.   Effect  of  Cavitation  on  Attenuator  Response 


t 

Values  of  Ap    /p 

Kmax  *s 

Ph     Ps 

GFC 

REL 

psia 

psia 

* 

N.C. 

W.C.      N.C.  \     W.C. 

30 

1000 

0.21 

0.69      0.71 

0.76 

280 

750 

; 

0.44 



0.45      0.82 

0.82 

rN.C 


no  cavitation 


**W.C. 


with  cavitation 


It  is  postulated  that  the  relation  between  gage  pressure 
p   and  thickness  L  is  p   =  C,  (L  /L-l). 
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Results  in  Table  1  show  that  cavitation  has  little  effect 
on  the  performance  of  the  REL  considered  at  either  hydrostatic 
pressure.   This  is  also  true  for  the  GFC  at  p,  =  280  psi.,  but 
there  is  severe  performance  degradation  at  p,  =  50  psia.   Note, 
however,  that  performance  remains  better  than  that  of  the  REL. 
If  the  hydrostatic  pressure  is  increased  significantly  above 
280  psia.,  maintaining  the  relation  p-,  +  p   =  const.,  cavitation 
will  not  occur. 

5.   Conclusions 

Both  the  displacement  formulation   and  the  displacement 
potential  formulation   have  been  shown  to  produce  acceptable 
results  when  applied  to  the  Bleich-Sandler  example.   It  is 
anticipated  that  either  formulation  will  provide  a  workable 
basis  for  solving  three-dimensional  axisymmetric  problems. 
The  scalar  displacement  potential  will  lead  to  a  much  smaller 
computer  storage  requirement  and  processing  time,  but  it  may 
not  be  easy  to  fit  it  to  the  framework  of  an  existing  program 
such  as  NASTRAN,  MARC,  or  NONSAP .   No  insurmountable  difficulty 
is  anticipated  in  using  the  displacement  formulation  with  one 
of  these  programs. 
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Appendix  A.   Governing  Equations 

Equations  are  derived  here  in  a  form  independent  of  the 
number  of  spatial  dimensions. 

Newton's  Second  Law:    p5_  =  -Vp  +  £.  (Al) 

In  the  above: 

p   =   fluid  density; 

5_  =   particle  displacement  vector; 

V   =   gradient  operator; 

p   =   fluid  absolute  pressure; 

£   =   body  force  per  unit  volume. 
Note  that  the  underline  is  used  to  denote  a  vector  quantity. 
Differentiation  with  respect  to  time  is  denoted  by  a  superior 
dot  and  the  convective  contributions  to  the  material  derivative 

are  neglected. 

2 
Bilinear  Constitutive  Law:   p   =   -c"pe,  e  <_   0; 

2  2  (A2) 

p   =   -6  c  pe ,  e  >  0 . 

Here : 

c   =   acoustic  velocity  in  fluid; 

e   =   dilatation. 

2 
Note  that  c  p  is  the  bulk  modulus  of  the  fluid.   For  the  bilinear 

fluid  model  8  is  chosen  as  positive  and  small  compared  with  unity 

The  limiting  condition  of  zero  pressure  in  the  cavitated  region 

corresponds  to  3  =  0 . 

Geometric  identity:   e   =   V*6_  ,  (A3) 

where  the  dot  denotes  the  scalar  product. 
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It  is  possible  to  choose  a  single  dependent  variable 
such  as  p  and,  through  suitably  chosen  manipulations,  demon- 
strate that  p  obeys  the  wave  equation  in  the  uncavitated  fluid 
and  a  modified  form  with  3c  in  place  of  c  in  the  cavitated 
region(s) .   Thus 

2  2 

p   =   cVp,   p^O; 

777  (A4J 

p   =   3  c  V  p,  p  <  0. 

A  more  enlightening  approach  which  focuses  attention  on 
the  sequential  steps  in  time  integration  of  the  governing 
equations  uses  auxiliary  dependent  variables  and  a  set  of 
three  equations.   For  this  purpose  we  first  define  some  addi- 
tional dependent  variables  and  then  summarize  four  separate 
formulations . 

Definitions .   In  our  applications  the  body  force  f  appearing 
in  (Al)  may  be  expressed  as 

f   =   Vph  ,  CA5] 

where  p,  is  the  hydrostatic  component  of  fluid  pressure. 

It  is  useful  to  introduce  two  density  weighted  variables: 

r   =   p5_  ,  (A6) 

s   =   -pe.  (A7) 

We  also  introduce  two  similarly  weighted  potential 
functions ; 

y    =    i,  (as) 

V^   =   r.  (A9) 

Henceforth  we  omit  explicit  reference  to  the  density  factor  and 
refer  to  r  as  displacement,  s  as  condensation  (Lamb's  usage, 
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see  Ref .  10)  ,  <£  as  velocity  potential  and  ip  as  displacement 
potential . 

Using  s,  the  bilinear  constitutive  law  is  rewritten  as 

2  n 

p   =   c  s ,      s    0 ; 

2  2  (A10) 

p   =   3  c  s ,    s<0. 

In  the  computational  stage  a  further  simplification  is  effected 
by  choosing  length  and  time  units  such  that  c  =  1. 

r  Formulation.   Using  (A5)  and  (A6) ,  (Al)  becomes 

r   =   -V(p-ph).  (All) 

Using  also  (A7) ,  (A3)  becomes 

s   =   -V-r.  (A12) 

When  applicable  initial  and  boundary  conditions  are  prescribed, 

the  r  formulation  allows  the  following  calculation  sequence: 

1.  Using  present  values  of  p  and  p,  ,  calculate  r  from  (All). 

2.  Using  a  suitable  time  integration  algorithm  and  the 
current  values  of  r  and  r,  find  new  values  of  r  and  r_  after 
one  time  step. 

3.  Use  (A12)  to  find  new  values  of  s. 

4.  Find  corresponding  new  values  of  p  from  (A10) . 

5.  Return  to  Step  1  with  new  values  of  p  and  repeat  the 
sequence  as  many  times  as  needed. 

p  Formulation.   Determining  the  divergence  of  both  sides  of 

(All),  then  calculating  the  second  time  derivative  of  each  side 

of  (A12)  and  substituting  in  the  preceding  result  gives 
2 


s 


VCp-Ph).  (A13) 
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It  appears  as  if  a  sequential  use  of  (A13)  and  (A10)  in  a 
fashion  paralleling  that  described  above  for  the  r  formulation 
would  allow  tracking  the  time  history  of  p.   The  process  is 
workable  in  the  absence  of  cavitation  and  has  been  success- 
fully applied  to  a  variety  of  problems  (e.g.,  see  Refs.  7-9). 
In  such  applications  the  variable   s   is  superfluous  and  the 
first  of  (A4)  suffices.   Reasons  for  the  failure  of  this  for- 
mulation in  a  cavitated  region  are  discussed  in  Art.  2.5. 


.=> 


<j)  Formulation.   Using  (A8)  and  (All)  we  may  deduce  the  result 
i      =   ph  -  p.  (A14) 

Further,  using  (A8)  and  (A12),  we  may  find 

s   =   -724>.  (A15) 

These  two  equations,  followed  by  (A10) ,  may  be  employed  sequen- 
tially and  repetitively  to  construct  a  time  marching  solution. 
Although  cf> ,  like  p,  satisfies  the  wave  equation  (A4)  ,  such  a 
reduction  of  the  equations  still  requires  calculation  of  s 
by  integrating  (A15)  to  distinguish  cavitated  regions. 

ip  Formulation.   Using  (A8)  and  (A9)  we  may  transform  (A14)  into 
J   =   ph  -  p,  (A16) 

and  transform  (A15)  into 

s   =    -V24>.  (A17) 

These  two  equations,  followed  by  (A10),  also  may  be  used 
sequentially  and  repetitively  to  find  the  time  history  of  p. 
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Appendix  B.   Initial  and  Boundary  Conditions 

Initial  Values  of  j>  and  j .   Considered  here  are  the  initial 
conditions  for  an  uncavitated  region  with  hydrostatic  pressure 
p-,  and  a  dynamic  pressure  p-   resulting  from  a  wave  travelling 
in  the  negative  x  direction.   Thus,  at  time  t: 


p(x,t)   =   ph(x)  +  pin(x,t) . 
Our  immediate  concern  is  with  conditions  at  t  =  0 


*, 


xx 


p3u/3x 


■p(x,0)/c*. 


(Bl) 
Now 
(B2) 


Integrating  twice  gives  the  result 

x   rX 


<Kx,0)   =   pxu(0,0)   -  -^ 

c 


0 


p(C,0)dCdX,     (B3; 


0 


where  the  choice  $(0,0)  =  0  is  arbitrary.   For  evaluation  of 
$(x,0)  we  begin  with 


¥, 


pu 


(B4) 


The  particle  velocity  is  induced  by  the  incoming  wave  and  is 
given  by 


Substituting  (B5)  into  (B4)  and  integrating 

x 


$(x,0)   =   pcu(0,0) 


Pin(C,0)dc 


(B6) 


0 


The  choice  of  $(0,0)  =  pcu(0,0)  is  useful  in  connection  with 
the  radiation  boundary  condition  considered  in  the  next  article 
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The  initial  conditions  given  by  (B3)  and  (B6)  are  based 
on  an  incoming  pressure  wave  in  uncavitated  fluid.   The  modi- 
fications required  to  deduce  initial  conditions  for  the  water- 
hammer  problem  are  obvious  and  are  not  detailed  here. 

Radiation  Boundary  Condition.   Representation  of  a  semi  -  infinite 
region  by  the  finite  element  method  requires  some  strategem  for 
truncating  the  discretized  region.   The  device  employed  here  is 
an  extension  of  the  radiation  boundary  condition  originally 
introduced  in  Ref.  7  and  successfully  employed  in  Ref s .  8  and  9. 
The  relations  used  are  based  on  the  d'Alembert  solution  to  the 
wave  equation.   Thus,  for  an  incoming  wave: 

*inC*,t)   =   fCx+ct).  (B7) 

Similarly,  for  an  outgoing  wave: 

*out(x,t)  =   gCx-ct).  (B8) 

For  our  problems  we  may  write 

\b  =       \b,     +  ib.       +    \b       ..  (B9) 

where  ip,  is  contributed  by  the  hydrostatic  pressure.   If  we 
choose  to  terminate  the  region  at  x  =  x   (the  radiation 
boundary),  we  require  \b ,  (x  ,t)  for  our  boundary  condition. 
Using  (B7)  ,  (B8) ,  and  (B9)  it  is  readily  established  that 

♦•x    =   *h,x  +  2*in,x  -  */c-  tBIO) 

By  rather  obvious  extensions  of  the  manipulations  leading  to 

(B3)  ,  the  needed  values  of  i|>,    and  fy .  may  be  found.   The 

n  j  a       in  ,  a 

value  of  ip    is  generated  in  the  solution  process. 
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Appendix  C.   Time  Integration  and  Artificial  Damping 

Time  Integration  Algorithm.   Prior  experience  with  transient 
wave  propagation  studies  by  the  finite  element  method  (e.g., 
see  Ref.  11)  established  the  desirability  of  using  a  time 
integration  algorithm  which  effectively  introduces  damping  that 
increases  with  modal  frequency.   Also  desirable  was  a  method 
explicitly  designed  for  second-order  equations.   Two  methods 
known  to  meet  these  requirements  are  the  Houbolt  method  (Ref.  12) 
and  the  Wilson  9  method  (Ref.  13) .   Both  of  these  algorithms 
can  be  unconditionally  stable.   The  Houbolt  method  introduces 
greater  spurious  damping  than  Wilson's  (Ref.  14).   This  advan- 
tage is  offset  by  the  fact  that  the  Houbolt  method  approximates 
the  second  time  derivative  by  fitting  a  cubic  polynomial  to 
four  equally  spaced  ordinates.   Because  of  the  discontinuities 
inherent  in  the  cavitation  problem  the  Wilson  method,  which 
utilizes  only  two  adjacent  ordinates  for  each  time  step,  was 
chosen. 

The  nonlinearity  of  the  governing  equations  in  the  neigh- 
borhood of  a  cavity  boundary  necessitates  a  nonstandard  appli- 
cation of  the  Wilson  method.   Using  8  =  1.4  the  method  assumes 
that  f  is  a  linear  function  of  time  from  the  current  instant 
for  a  duration  1.4h,  where  h  is  the  time  step.   For  this  appli- 
cation, an  initial  estimate  of  the  forward  value  of  i>   was  based 
on  linear  extrapolation.   The  estimate  was  improved  by  iteration 
before  proceeding  to  the  following  time  step.   An  effect  of 
using  this  strategem  was  to  introduce  a  limit   on  the  maximum 
usable  time  step  (i.e.,  to  sacrifice  unconditional  stability). 
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Artificial  Damping.  Initial  solutions  using  the  Wilson  method 
showed  both  temporal  and  spatial  oscillations  of  pressure  fol- 
lowing passage  of  the  shock  front  in  the  water-hammer  problem. 
Since  no  such  behavior  is  shown  by  the  exact  solution,  it  is 
clearly  a  numerical  artifact.  To  suppress  the  unwanted  oscil- 
lation, damping  was  introduced  into  the  governing  equations. 
The  mechanism  chosen  was  to  modify  (A16)  to  read 

'i      =  Ph  "  P  "  ns.  (CI) 

The  coefficient  n  appearing  in  (CI)  was  chosen  by  cut-and-try. 
The  needed  value  of  s  is  calculated  from  (A17)  by  differentiat- 
ing with  respect  to  time. 
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